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Erosion by ion-beam sputtering (IBS) of amorphous targets at off-normal incidence frequently 
produces a (nanometric) rippled surface pattern, strongly resembling macroscopic ripples on aeolian 
sand dunes. Suitable generalization of continuum descriptions of the latter allows us to describe 
theoretically for the first time the main nonlinear features of ripple dynamics by IBS, namely, 
wavelength coarsening and non-uniform translation velocity, that agree with similar results in ex- 
periments and discrete models. These properties are seen to be the anisotropic counterparts of 
in-plane ordering and (interrupted) pattern coarsening in IBS experiments on rotating substrates 
and at normal incidence. 
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Ever since their earliest observation [1], the pro- 
duction of ripples on the surfaces of amorphous tar- 
gets subject to ion-beam sputtering (IBS) at inter- 
mediate energies, has been fascinating due to the 
similarities with macroscopic ripples, like those pro- 
duced underwater [2], or on the surface of aeo- 
lian sand dunes [3]. Beyond the morphological re- 
semblance, IBS ripples share many other proper- 
ties with e.g. aeolian ripples, such as wavelength 
coarsening and pattern translation with time [4, 5]. 
Remarkably, while typical wavelengths of the lat- 
ter are above 1 cm, the periodicity of IBS ripples 
is in the 100 nm range [6], these patterns having 
gained increased interest recently for applications in 
Nanotcchnology, ranging from optoelectronic to cat- 
alytic [7]. IBS ripples arc produced on a wide class 
of substrates, from amorphous or amorphizable (sil- 
ica, Si, GaAs, InP) to metallic targets (Cu, Au, Ag) 
[6] . In view moreover of their implied loss of in-plane 
symmetry, IBS ripples provide interesting instances 
of systems hosting a competition between pattern 
forming and disordering mechanisms [8] . 

A successful description of the main features of 
IBS ripples was provided by Bradley and Harper 
(BH) [9], based on Sigmund's linear cascade ap- 
proximation of sputtering processes in amorphous 
or polycrystalline targets [10]. The linear equation 
derived by BH describes satisfactorily some proper- 
ties of IBS ripples, such as their alignment with the 
ion beam as a function of the incidence angle to tar- 
get normal 9 [wave vector parallel (perpendicular) to 
the projection of the ion beam for 6 < 8 C (9 > 6 C ), 
for some threshold 9 C \. Other features, such as ripple 
stabilization or wavelength dependence with ion en- 
ergy or flux, required non-linear extensions of BH's 
approach [11, 12], leading to an anisotropic gener- 



alization of the well known Kuramoto-Sivashinsky 
(KS) equation [8, 11]. However, a notable limitation 
of the anisotropic KS (AKS) equation is its inabil- 
ity to predict ripples that coarsen with time, con- 
tradicting observations in many experiments and/or 
discrete models of IBS (see [5, 13] and Refs. therein). 

In this Letter, we introduce a "hydrodynamic" 
model [3, 14, 15] for IBS ripple production at off- 
normal incidence. Time scale separation between 
microscopic processes (collision cascades, surface dif- 
fusion) and pattern evolution allows us to derive an 
improved equation for the surface height. The new 
non-linear terms appearing allow for ripple coars- 
ening and pattern translation with non-uniform ve- 
locity, as seen in experiments and discrete models. 
Our theory has both the AKS equation [11] and the 
normal incidence hydrodynamic theory [15, 16] as 
particular limits, and enables analysis of the impor- 
tant case of rotating substrates [17]. In addition, our 
model may be important also to the context of rip- 
ples on aeolian sand dunes, where the standard ID 
approximation requires validation, fully anisotropic 
2D models being scarce [3] , as incidentally occurs in 
many other contexts within Pattern Formation [18]. 

During IBS of amorphous or semiconductor sub- 
strates, in which the ions amorphizc the subsurface 
layer, incident ions lose their energy through ran- 
dom collision cascades in the bulk [10]. Near-surface 
atoms receiving enough energy and momentum to 
break their bonds are in principle sputtered away, al- 
though they may join the current of surface adatoms 
that are available to other relaxation mechanisms, 
such as surface diffusion, before incorporating back 
to the solid bulk. Within the so-called "hydrody- 
namic" approach to aeolian sand dunes [3] and ion- 
sputtered surfaces [14, 15], we consider two coupled 
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fields, namely, the thickness of the mobile surface 
adatoms layer, i?(x, t), (related with the density of 
mobile adatoms through the atomic volume) and the 
height of the bombarded surface above a reference 
plane, /i(x, t). Their time evolution is provided by 

d t R = (i - &r ex - r ad + dv 2 r, (i) 

d t h = -T ex + T ad , (2) 

where T ex and T a d are, respectively, rates of atom 
excavation from and addition to the immobile bulk, 
(1 — 4>) = <f> measures the fraction of eroded atoms 
that become mobile, and the third term in Eq. (1) 
describes motion of mobile atoms along the surface 
as due to isotropic thermal diffusion (D being a con- 
stant for amorphous materials). Even if all eroded 
atoms are sputtered away (<p = 1), we assume a non- 
zero average fraction of mobile atoms, R eq - 

Under the assumption that nucleation events are 
more likely in surface protrusions, in analogy to the 
Gibbs-Thompson relation we have 

Tad = 7o [R[l + l2 X d 2 x h + l2y d 2 y h) - R eq ] , (3) 

where 70 is the mean nucleation rate for a flat surface 
(on the xy plane) and 72^, 72^ > describe variation 
of the nucleation rate with surface curvatures. Note 
that, in (3), the full thickness of the mobile adatoms 
layer is affected by the shape of the surface. 

The rate at which material is sputtered from the 
bulk depends on the angle of incidence, ion and sub- 
strate species, ion flux, energy, and many other ex- 
perimental parameters. If the beam direction is in 
the xz plane, we have, following [11, 15] 

r e:E = a [l + a lx d x h + ct 2x d 2 x h + a 2y d 2 h 
+a 3x {d x h) 2 +a 3y (d y h) 2 -(d x h){a 4x dlh+a 4y d 2 h)], 

(4) 

where parameters reflect the dependence of T ex on 
the local shape of the surface [19], as described by 
more microscopic derivations such as BH or general- 
izations thereof [11, 20]. Analogous, but not equal, 
geometrical couplings to the driving occur in aeo- 
lian sand dunes [3], or in growth onto amorphous 
substrates [21]. Note the loss of reflection symmetry 
in the x direction, but not in the y direction. For 
a planar surface, atoms are sputtered from the bulk 
in a typical time of order ajT 1 . 

The main difference between former models [9, 11, 
12] and our present model, Eqs. (l)-(4), is that, in 
the latter, eroded material is allowed to redeposit 
locally, and there is an implicit viscous flow [22] in 
the amorphized layer through the evolution of R. 
These additional mechanisms are seen below to in- 
duce richer pattern dynamics than in [9, 11, 12]. 



The linearized Eqs. (l)-(4) have solutions R l = 
i?'e lkx+Wkt , h l = /i'e lkx+ " kt , with a dispersion re- 
lation, tJk; given in the long wavelength limit by [19] 

T^e(wk) = e(/)~fo(a 2 xkl + ct 2y k 2 ) - e 2 cj>^ a 2 lx k 2 x 

-R eq D(k 2 x + kl){ l2x k 2 x + l2y k 2 y ) + <D{ek% (5) 
Xm(w k ) = -e(jyy a lx k x + £>(fc 3 ), (6) 

where e = ao/70 is a dimensionless parameter; the 
erosion rate being much smaller than the nucleation 
rate, e <C 1 for typical experiments. Eq. (6) is a sim- 
ple consequence of the asymmetry in the x direction, 
induced by the incoming flux. 

The surface morphology is dominated by the peri- 
odic pattern with wave vector ]t max making Eq. (5) 
a positive maximum. It can be shown [19] that ]f max 
is oriented along the x or y directions, as observed 
experimentally [6]. Close to the instability thresh- 
old, before nonlinear terms are no longer negligible, 
one has k™y X ~ e 1 / 2 . Substituting this into Eqs. (5) 
and (6) provides us with estimations of the typical 
time and length scales of the pattern, that we em- 
ploy to rescale X = e 1/2 x, Y = e 1/2 y, T 2 = e 2 t and 
Ti = e 3 / 2 £, and perform a multiple scale expansion 
of Eqs. (1), (2), in which R can be adiabatically elim- 
inated. To lowest non-linear order, we get [19, 23] 

d t h = ~fd x h + {-"i 9fh + Aj 1 ) (d<h) 2 + ^d 2 d x h 

i=x,y 

+ ^ dxh )(dfh)}+ J2 {-^ijdfd'h + X^dUdjh) 2 }; 

hj=x,y 

7 = -e07 O ai K , v x = ecj)-i a a 2x - e 2 (j)(jry a 2 lx , 
v y = e4>j a 2y , \r' = -e(jry a 3i , 
Cli = e(4>D - (j)R e qlol2i)cti Xl & = e07oa 4i , 
Kij = DR eq ~f 2l + e[Dcf>(~f 2l - a 2 j) + HoRe q J 2 ia 2j ] , 
^ij ) = 6 ~ 4>Re q l0l2i] a 3j ■ (7) 

Eq. (7) generalizes the AKS type equations obtained 
[9, 11, 12] within BH approach to IBS. While shar- 
ing the same reflection properties in the x and y 
directions and most of the terms on the rhs, both 

equations differ crucially by the presence here of the 
( 2 ) 

\\j nonlinearities. Moreover, in the absence of rede- 
position (<j) = 1), A^ 1 "* and A^ have the same signs, 
making Eq. (7) nonlincarly unstable, as in the BH 
case [15, 24, 25]. Note, the linear dispersion rela- 
tion of (7) matches Eqs. (5) and (6) above. Un- 
der normal incidence, parameters are isotropic and 
o-\x = ctix = ciiy = 0, Eq. (7) reducing to that ob- 
tained in [15, 24], and in [21]. As a final remark, 
let us quote the form of Eq. (7) for sample rota- 
tion around the z axis during bombardment (see e.g. 
[17, 26]). Dynamics of h are given by a different 



isotropic limit of Eq. (7), namely, 

dth = -v,V 2 h - JC^h + X^iVh) 2 + X (2 W 2 (Vh) 2 
+ \[ 3) V- [(V 2 /i)V/i]] ; (8) 
Ur = {u a + u y )/2, A« = (A«+A«)/2, 
A (2) _ I V A (2) A^ = i V A (2) <5- ■ - X {2) 

/C r = (3/C T ^ + 3/Cy y ~\~ K- X y ~\~ K-y x) / 8 , 

with parameters (susbscript "r" denotes "rotating"). 

To the best of our knowledge, Eq. (7) is new, and 
indeed has a rich parameter space. Numerical in- 
tegration (not shown) within linear regime retrieves 
all features of the ripple structure as predicted by 
the BH theory, i.e., dependence of the ripple wave- 
length with linear terms, and ripple orientation as 
a function of 8. Entering the nonlinear regime, and 
as occurs in the AKS equation and its generaliza- 
tions [11, 12], nonlinearities lead to saturation 
of the pattern with constant wavelength and ampli- 
tude. In absence of these terms, the ripple wave- 
length grows indefinitely as £(t) ~ i" with n = 1/2 
[27, 28] until a single ripple remains in a finite sim- 
ulation domain. As an anisotropic generalization of 
the ordering process observed for normal incidence 
[15, 16], pattern coarsening requires the presence 

(2) 

of A),-, whose magnitude and mathematically cor- 
rect sign [25] are due to describing redeposition by 
means of the additional field R. When the values of 
these coefficients increase relative to A^ , coarsen- 
ing stops later, and the amplitude and wavelength 
of the pattern also increase. The coarsening expo- 
nent n will take an effective value that will be larger 
the later coarsening stops, and may depend on sim- 
ulation parameter values. For instance, we show in 
Fig. 1 snapshots of a numerical integration of Eq. (7) 

(2) 

for relatively large A^ ■ . The apparent coarsening is 
quantified in the plot of £{t) shown in Fig. 2(a), com- 
patible (after transient effects analogous to those in 
[27]) with n = 0.19. Note the saturation of ripple 
wavelength at long times, together with saturation 
in amplitude, as shown in the plot of the surface 
roughness (rms width) W(t) in Fig. 2(6). These re- 
sults are similar to those obtained experimentally 
for IBS of Si in [4]. For a different experiment on Si, 
precise measurements of the coarsening exponent [5] 
yield n = 0.50(4), no saturation having been ob- 
served in this system, as we would expect. Here, 
the dispersion velocity of the pattern was also mea- 
sured, finding that it decays with ripple wavelength 
or, equivalently, with time. We have also observed 
the same trend in the dispersion velocity of the rip- 
ples shown in Fig. 1. Experiments also exist, e.g. 
for IBS of Si, in which ripple coarsening is absent or 
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FIG. 1: Numerical integration of Eq. (7) using u x — 1, 

u y = o.i, n x =i,n y = 0.5, & = o.i, a! 1 ' = o.i, a^ = 

5, Afj = -5, K,i,j = 1. Top views for t = 10, 106, 953 
(a), (b), and (c), respectively. Insets are height autocor- 
relations. All units are arbitrary, (d): Top to bottom, 
side cuts of (c), (fe), (a), and grooves obtained for sign- 
changed X^c) v and X[ 2 j. Curves have an artificial offset. 
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FIG. 2: Ripple wavelength £(t) and surface roughness 
W(t) for parameters as in Fig. 1. The straight line rep- 
resents £ ~ t ' 19 . All units are arbitrary. 

(2) 

residual [29], that would correspond to smaller A^- 
values in (7), see e.g. an example in Fig. 3, where 
£(t) ~ log£ approximately. Ripple coarsening has 
been also observed in a Monte Carlo model of IBS 
[13], in which rules implement Sigmund's theory. To 
our purpose, the main conclusion of this study is 
the correlation between an increasing £(t) and non- 
uniform dispersion velocity, and on the parameter- 
dependent values of the coarsening exponent. Let us 
note that, on the experimental side, values of n show 
a large scatter in the literature (see Refs. in [6, 13]). 
Additional non-linear effects can be described by Eq. 
(7). A first one is production of grooves (as opposed 
to ripples), due to the loss of up-down symmetry in- 
duced by quadratic nonlinearities. Indeed, by chang- 
ing the sign of A^ 1 , grooves replace ripples, see Fig. 
1(d); this calls for systematic experimental explo- 
ration. A second effect is related with cancellation 
modes (CM), known in the AKS equation [12, 30] 
and other models of IBS [14, 15, 24]. These are lin- 
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FIG. 3: Numerical integration of Eq. (7) for parameters 
as in Fig. 1, except for = 1. Top views for t = 
10,106,953 (a), (b), and (c), respectively, (d): Surface 
roughness W(t) for the same simulation. Inset: £(t) for 
the same system, showing marginal coarsening. All units 
are arbitrary. 



early unstable modes for which nonlinearities cancel 
one another exactly. While for equations like the 
one considered in [15, 21] CM affect well-posedness, 
anisotropic systems [12, 30] may remain better de- 
fined in the presence of CM. In our case, if reflection 
symmetry breaking terms are neglected in Eq. (7) 
and BH parameters are used, the same CM occur as 
in the AKS equation. Additional CM ensue between 
the and terms for appropriate relative signs 
[21, 25]. We have verified numerically that the full 
Eq. (7) breaks down for the latter CM, but can sup- 
port AKS-type CM. That these solutions are phys- 
ically realizable or are artifacts of the small slope 
approximation made, remains to be assessed. Use- 
ful information on this issue might come from field 
experiments in aeolian sand dunes. 

In closing, we mention IBS of metals as an im- 
mediate experimental domain to which the above 
results may be relevant [6], albeit differing in the 
degree of universality. There, however, the correct 
extension of BH theory is not yet clear, nor is its 
importance relative to anisotropic surface diffusion. 
We have taken preliminar steps in this context [20] , 
and expect to make progress in this direction soon. 
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